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Abstract. We investigate the fitness advantage associated with the robustness of a 
phenotype against deleterious mutations using deterministic mutation-selection models 
of quasispecies type equipped with a mesa shaped fitness landscape. We obtain analytic 
results for the robustness effect which become exact in the limit of infinite sequence 
length. Thereby, we are able to clarify a seeming contradiction between recent rigorous 
, work and an earlier heuristic treatment based on a mapping to a Schrodingcr equation. 

Qpi We exploit the quantum mechanical analogy to calculate a correction term for finite 

sequence lengths and verify our analytic results by numerical studies. In addition, we 
investigate the occurrence of an error threshold for a general class of epistatic landscape 
^ ' and show that diminishing epistasis is a necessary but not sufficient condition for error 

C*") , threshold behavior. 

q 

O ■ 1. Introduction 

<3\ ' 

In current evolutionary theory, the concept of robustness, referring to the invariance 
of the phenotype under pertubations, is of central importance [H [2] . Here we address 
specifically mutational robustness, which we take to imply the stability of some biological 
function with respect to mutations away from the optimal genotype. To be precise, 
suppose the genotype is encoded by a sequence of length L, and the number of 
mismatches with respect to the optimal genotype is denoted by k. Robustness is then 
quantified by the maximum number of mismatches ko, that can be tolerated before the 
fitness of the individual falls significantly below that of the optimal genotype at k = 0. 
This situation arises e.g. in the evolution of regulatory motifs, where the fitness is a 
function of the binding affinity to the regulatory protein [31 H] . Assuming that the fitness 
is independent of k both for k < k and for k > k , the fitness landscape has the shape 
of a mesa parametrized by its width ko and height wq [5]. 

We consider deterministic mutation-selection models of quasispecies type, which 
describe the dynamics of large (effectively infinite) populations [6J. We analyse the 
stationary state of mutation-selection balance, focusing on the dependence of the 
population fitness on the parameters ko and wq. This allows us to identify the conditions 
under which a broad fitness peak of relatively low selective advantage outcompetes a 
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higher but narrower peak [7J, a phenomenon that has been referred to as the survival 
of the flattest [HI [9J [TQl |TT] [12]. Our central aim is to obtain analytic results for the 
robustness effect that become exact in the limit of long sequences. In particular, we 
want to clarify whether the selective advantage is a function primarily of the relative 
number of tolerable mismatches xq = ko/L, or of the total number of mismatches ko- 

Robustness in the sense described above is a special case of epistasis, which refers 
more generally to any nonlinear relationship between the number of mutations away 
from the optimal genotype and the corresponding fitness effect [T3] . A simple way to 
parametrize epistasis is to let the loss of fitness vary with the number of mismatches 
as k a , such that the non-epistatic case a = 1 separates regimes of synergistic (a > 1) 
and diminishing (a < 1) epistasis [HI [15]. An important problem in previous work 
on mutation-selection models has been to identify the conditions under which epistatic 
fitness landscapes display an error threshold, a term that refers to the discontinuous 
derealization of the population from the vicinity of the fitness peak as the mutation 
rate is increased beyond a critical value [Bl HE]. Improving on earlier work that found 
that only landscapes with diminishing epistasis (a < 1) have an error threshold, we 
derive here the more stringent condition a < 1/2 on the epistasis exponent. 

1.1. Organization of the article 

We base our work on two complementary analytic approaches. First, recent progress 
in the theory of mutation-selection models j5j [TBI El [13 HSl EDI ED [22] provides 
an expression for the population fitness in terms of a maximum principle (MP) that 
becomes exact when the limit L — > oo is performed keeping the ratio xq = ko/L fixed. 
Second, Gerland and Hwa (GH) [3] have used a drift-diffusion approximation to map 
the mutation-selection problem onto a one-dimensional Schrodinger equation which is 
then analyzed with standard techniques. 

Our work was initially motivated by the observation of a discrepancy between the 
two approaches: Whereas the MP predicts that the selective advantage of a broad mesa 
should vanish when the limit L — > oo is taken at fixed ko, in the GH approach a finite 
selective advantage is retained in this limit, which depends on the absolute value of ko 
rather than on xq. After introducing the model and briefly reviewing the results of the 
MP approach in section [2], we therefore provide a detailed discussion of the drift-diffusion 
approximation used by GH in section [3j We emphasize that it amounts to a harmonic 
approximation, and show how it can be improved in such a way that the results based 
on the MP are recovered. 

The mapping to one-dimensional quantum mechanics is nevertheless useful, as it 
allows us to derive the leading finite size correction to the population fitness. As a 
consequence we find excellent agreement between the analytic predictions and numerical 
solutions of the discrete mutation-selection equations. In section H] we consider the 
selection transition in a two-peak landscape first studied by Schuster and Swetina 
[71 [23], in which the population shifts from a high, narrow fitness maximum to a lower 
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but broader peak with increasing mutation rate. The occurrence of this transition 
is an indicator for the superiority of robustness over fitness in certain parameter 
regimes. In section [5] we use the MP approach to derive the critical value of the 
epistasis exponent a and verify our prediction by numerical calculations. Finally, 
some conclusions are presented in section El Details of the derivation of the improved 
continuum approximation and the generalization to arbitrary alphabet size can be found 
in two appendices. 



2. Mutation-selection models and the maximum principle 

We consider the simplest case of binary sequences and adopt continuous time dynamics 
of the Crow-Kimura type, in which the mutation and selection terms act in parallel 
[Sj. Point mutations occur at rate /i, and the (Malthusian) fitness is assumed from the 
outset to be a function w k only of the Hamming distance k to the optimal sequence at 
k = 0. The population structure is described by the fraction P k (t) of individuals with 
k mismatches, which satisfies the evolution equation 
dP 

-± = (w k - w)P k + fx(k + l)P k+1 + n(L-k + l)P fc _! - fiLP k . (1) 
dt 

with 1 < k < L — 1 and obvious modifications for k = and k = L. The nonlinearity 
introduced by the mean fitness w(t) = J2 k w k P k can be eliminated by passing to 
unnormalized population variables [6j [IM]. At long times the population distribution 
therefore approaches the principal eigenvector P£ of the linear dynamics, which is the 
solution of the eigenvalue problem 

KP* k = (w k - vL)P* + fj,(k + 1)P; +1 + fi(L-k+ l)P fc *_! (2) 

with the maximal eigenvalue A. This eigenvalue is equal to the long-time limit of 
the mean population fitness w, and it is the main quantity of interest in this paper. 
Depending on the context we will refer to A as the mean population fitness, the 
population growth rate, the principal eigenvalue of the mutation-selection matrix defined 
by (EJ) or the ground state energy of the corresponding quantum mechanical problem, 
to be defined in subsection 13.21 

A considerable body of work has been devoted to the solution of (|2]) for large L. In 
order to obtain nontrivial behavior in the limit L —>■ oo, it is necessary to either scale 
the mutation rate ~ 1/L or the fitness ~ L. We adopt here the first choice and take 
L — > oo, n — ► with 7 = \iL fixed. If, in addition, the fitness landscape w k is assumed 
to depend only on the relative number of mismatches, such that 

w k = f{x), x = k/L (3) 

the principal eigenvalue in (TSJ) is given, for L —>■ oo, by the solution of a one-dimensional 
variational problem as pH QH Mi UM EH M\ 

A = max {/ 0) - 7 [1 - 2^(1 -£)]}; (4) 
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see subsection 13.51 and Appendix A for a heuristic derivation, and Appendix B for the 
generalization to arbitrary alphabet size. Moreover, if f(x) is different iable the leading 
order correction to (ED takes the form 1 191 [20] 



AA = 2- =[1-^1-2/^)^-^)3/77], (5) 

z±j\i x c x c 

where x c is the value at which the maximum in (jlj) is attained. 

In the first part of this paper we focus on mesa landscapes of the form 

f w > : < k < k 
1 : k > k , 

where Wq is the selective advantage of the functional phenotype and k denotes the 
number of tolerable mismatches. Within the class of scaling landscapes ([3]), this is 
realized by setting 

f(x) = w o 0(x-x o ), (7) 

where 6 is the Heaviside step function and xq = k /L. Provided xq < 1/2, application 
of the maximum principle (j4j) yields 



A 



W -7(1 - 2^X0(1 -X )) : W > W C Q = 7(1 - 2^X0(1 - X )) 

: w < Wq. 



The value Wq of the selective advantage marks the location of the error threshold at which 
the population delocalizes from the fitness peak and the location x c of the maximum in 
(jl]) jumps from x c = x to x c = 1/2. 

The expression ([5]) is clearly not applicable to the discontinuous mesa landscape 
((?]). In fact we will show below that the leading order correction AA is of order L~ 2 / 3 
or L -1 / 2 rather than L _1 in this case. 



3. Continuum limit and the drift-diffusion equation 

3.1. Derivation and status 

A natural approach to analyzing ([T]) and (T5]) for large L is to perform a continuum limit 
in the index k. To this end we introduce e = 1/L as small parameter and replace the 
population variable Pk by a function 

(f>(x) = lim P xL . (9) 

L— >oo 

The fitness is taken to be of the general form ([3]). Expanding the finite differences in 
(J2]) to second order in e then yields the stationary drift- diffusion equation 

/0- e 7^(l-2x)0 + ^^ = A0. (10) 

This is identical to the equation obtained by GH [3], who however write it in terms of 
the unsealed variable k = Lx. 



Robustness and epistasis in mutation- selection models 



5 



Before proceeding with the analysis of fllOl) . some remarks concerning the accuracy 
of the second order expansion are appropriate. In the absence of selection (/ = 0) the 
principal eigenvalue in ( flOl) is readily seen to be A = 0, and the corresponding (right) 
eigenfunction is a Gaussian centered at x — 1/2, 

0o(x) ~exp[-(l-2x) 2 /2e]. (11) 

This is just the central limit approximation to the binomial distribution 

Pl = 2- L ([) (12) 



k / 

which solves (j2J) for w k = and A = 0. It is well known that the central limit 
approximation of ( fi"2l) is accurate in a region of size VX around k = L/2, but becomes 
imprecise for deviations of order L. An improved approximation is provided by the 
theory of large deviations [25], in which the ansatz 

P k ~ exp[-Lu(x)\ (13) 

is made to obtain an expression for the large deviation function u(x). In the context of 
mutation- select ion models, this approach has recently been introduced by Saakian [20J, 
who showed that it allows to derive the exact relation in a relatively straightforward 
manner (see Appendix A). Equivalent results can be obtained by continuing the 
expansion in (1101) to all orders in e and treating the resulting equation in a WKB-type 
approximation, which essentially corresponds to the ansatz f fl3|) . see [22] . 

We conclude that the drift-diffusion approximation ( flO|) can be expected to be 
quantitatively accurate only near the center x = 1/2 of the sequence space. We will 
nevertheless adhere to this approximation in following three subsections, because it 
allows us to make contact with the work of GH and to formulate the eigenvalue problem 
([2]) in the familiar language of one-dimensional quantum mechanics. In subsection 13.51 
we then show how to go beyond the second order approximation. 

3.2. Mapping to a one- dimensional quantum problem 

The key step in reducing ( TlOl) to standard form is to symmetrize the linear operator on 
the left hand side, thus eliminating the first-order drift term. This can be achieved by 
the transformation 

4>(x) = v^)^(x), (14) 

with 4>o(x) from (ITT]) , which leads to the stationary Schrodinger equation 

2 j2 

- € ~y^ + V ( X W = -( A - ^ ( 15 ) 

with the effective potential 

nr) = |(l-2x) 2 -/(*). (16) 

The latter consists of the superposition of a harmonic oscillator centered around x = 1/2 
with the (negative) fitness landscape. As pointed out in [5], the inverse sequence length 
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e plays the role of Planck's constant H, which implies that the case of interest is the 
semiclassical limit of the quantum mechanical problem. In particular, for e — > the 
ground state energy —A becomes equal to the minimum of the effective potential. We 
thus arrive at the variational principle 

A= max[/(x)-^(l-2x) 2 ], (17) 

£€[0,1] Z 

which is precisely the harmonic approximation (in the sense of a quadratic expansion 
around x = 1/2) of the exact relation (jlj). In this perspective the error threshold 
corresponds to a shift between different local minima of V(x), which become degenerate 
at the transition point. The transition is generally of first order, in the sense that 
the location x c of the global minimum jumps discontinuously. Within the harmonic 
approximation the transition for the mesa landscape occurs at 



Wo c = -(l-2xo)^^l--^j (18) 
when xo = k /L <C 1 . 

3.3. Semiclassical finite size corrections 

For small but finite e, quantum corrections to the classical limit (ITT)) have to be taken 
into account. If f(x) is smooth, the ground state wave function is localized near the 
minimum x c of the effective potential, and the shift in the ground state energy can be 
computed by replacing V(x) by a harmonic well, 

V(x) « V(x c ) + \v"(x c ){x - x c f = V(x c ) + ~[4 7 - f"(x c )}(x - a; c ) 2 .(19) 

Identifying I/7 with the mass m of the quantum particle [compare to (IT5)) ]. we see that 
this corresponds to a harmonic oscillator of frequency u = 2 7 a/1 — /"(o; c )/4 7 . The 
ground state energy eu/2, together with the shift e 7 on the right hand side of (IT5I) . thus 
gives rise to the leading order correction 

AA = 1[1- Vl-/"(x c )/4 7 ], (20) 

which coincides with (JSJ) evaluated for x c ~ 1/2. Similarly, the width of the wave 
function is given by[j] 



£=J^e72^= Vn — . (21) 

[8y/l-nx c )/±y]V* 

In the case of the mesa landscape ([T]), the potential near x c = x consists of a linear 
ramp of slope 

- a = V'(x ) = 2 7 (2x - 1) < (22) 

followed by a jump of height wq. For small e, the jump can be considered as effectively 
infinite (as the kinetic energy of the particle is then very small), and the corresponding 



I Note that, because of the factor V0o in (|14|) - this is not equal to the width of the stationary population 
distribution. 
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quantum mechanical ground state problem is standard textbook material [26J. ^From 
the solution we obtain the prediction 

AA = Zl {h 2 /2mf'^ = 2 1 / 3 z l7 (l - 2x ) 2/3 L- 2 / 3 + 0{L~ l ), (23) 

where Z\ ~ —2.33811... is the first zero of the Airy function. The scaling AA ~ L~ 2 / 3 
was already noted in [5] . The width of the wave function can be estimated to be of the 
order 

i ~ (ft 2 /™*) 1 / 3 ~ e 2 / 3 (24) 

in this case. 



3.4- The quantum confinement regime 

We are now prepared to make contact with the approach of GH [3]. Assuming from 
the outset that the maximal number of mismatches is small compared to the sequence 
length, 1 <C k <C L, they neglect the contribution 2x in the drift term on the left hand 
side of fflOl) . The linear operator can then be symmetrized by the transformation 

= e x/e *fj{x), (25) 

which is obtained from (fl4l by neglecting the terms quadratic in x in O - This leads to 
a Schrodinger problem similar to (fl5l) . but with a potential that differs from —f(x) only 
by the constant term 7/2. The error threshold is determined by the point at which the 
decay of the wave function ip(x) matches the exponential factor e x l e in (T25l) . such that 
4>{x) ceases to be normalizabkj§. For k$ 1 the location of the transition is found by 
GH to be 

<=i( i+ S)- ( 26) 

which depends on the absolute number of mismatches k , but is independent of L. 

To reconcile this with the result (|18|) . we note that the semiclassical approximation 
must break down when the width of the semiclassical wave function, as estimated in 
subsection 13.31 becomes comparable to the width of the potential well provided by the 
fitness function. For the discontinuous mesa landscape this occurs when 

£ ~ e 2 /3 „ XQ = eko A; ~ e -1 / 3 = L 1/3 . (27) 

For a mesa that is shorter than L 1//3 , the energy of the wave function is determined 
by its confinement on the scale Xq, and it can be estimated from standard quantum 
mechanical considerations to be of the order of h 2 /{mx1) ~ r je 2 /x 2 ] ~ l/k 2 .- For 
k L 1 / 3 this supersedes the contribution ~ k /L on the right hand side of ([TBI . We 
conclude, therefore, that the leading "quantum" correction to the " classical" eigenvalue 
A = wo — 7/2 is a negative contribution proportional to 7/&0, which leads to a 

§ This requires to decay on a scale of order unity in unsealed coordinates at the transition, which is 
actually inconsistent with the assumption of slow variation on the scale of the sequence index k that 
underlies the continuum approximation. 
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corresponding positive shift in Wq, in qualitative agreement with (I26p . For smooth fitness 
landscapes the breakdown of the semiclassical regime occurs already at k ~ L 1 / 2 , but 
the condition for the confinement energy contribution 7/ &o to dominate the ko / L-term 
in (HHD still reads k < L 1 / 3 . 



3.5. Beyond the harmonic approximation 

So far, we have worked in the harmonic approximation around x = 1/2, which breaks 
down near the boundaries x = and x — 1. However, to access the regime 1 -C k <C L 
considered by GH, an accurate treatment of the region of small 1 C 1 is clearly 
necessary. In Appendix A we show how the quantum mechanical treatment can be 
extended such that it become quantitatively valid over the whole interval < x < 1. 
Based on the considerations of [201 . we arrive at the modified Schrodinger equation 



d 



2 



- e 2 -f^x(l-x)—^ + [7(1 - 2y/x(l-x)) - /(x)J ^ = -AV>, (28) 
which differs from ( |T5|) in two respects. First, the potential (|T6l) is replaced by 



V m (x) = 7(1 - 2y/x(l-x)) - f(x). (29) 

In the asymptotic limit e — > the principal eigenvalue is obtained by minimizing Vf u u, 
which exactly recovers the maximum principle (j3J). Second, the mass of the quantum 
particle described by (j2SJ) becomes position dependent, 



m{x) = ^27^/a;(l - x)J , (30) 

which replaces the simple identification m = I/7 in the harmonic case. Inserting ( |29l ) 
and (130|) into the expression (j23|) for the finite size correction yields 

AA = 2- 1/3 z ie - 2/3 7 (l - 2£ ) 2/3 Ml - x )}- 1/6 . (31) 

For fixed xq this still scales as e 2 / 3 = L -2 / 3 , but when taking L — > 00 at fixed ko, such 
that £0 — > 0, we find instead that 

AA ^ 2- 1 / 3 2l7 x - 1/6 e 2 / 3 = 2- 1 / 3 ^ l7 A;- 1/6 L- 1 / 2 . (32) 

We next revisit the considerations of subsection 13.41 The width of the ground state 
wave function is of order £ ~ (h 2 /ma) 1 / 3 , where both m and a now diverge as Xq 1 ^ 2 for 
xo — > 0. Consequently (|24l) is replaced by 

£ ~ (6 2 x ) 1/3 = e^ /3 , (33) 

and we see that the condition £ ^> ek for the breakdown of the semiclassical 
approximation is never be satisfied. We conclude that the quantum confinement regime 
discussed in subsection 13.41 in fact does not exist, and hence the improved semiclassical 
expression ( !3TI) for the finite size correction is expected to remain valid for all ko and 
all L, provided that ko, L 3> 1. 
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Figure 1. Growth rate A as a function of the plateau width ko for two values of the 
plateau height wo = 0.5,0.95. The sequence length is L = 100 and the mutation rate 
per sequence is 7 = 1. The solution of the maximum principle together with the L -1 / 2 - 
correction term (including the position-dependent mass) provides the best agreement 
with the numerics. The numerical values of the growth rate have been obtained by 
(numerical) calculation of the largest eigenvalue of the matrix defined by equation 




Figure 2. Critical plateau height as a function of the plateau width ko. The sequence 
length is L = 500 the mutation rate per sequence 7 = 1. The solution of the maximum 
principle together with the L _1 / 2 -correction term provides the best agreement with the 
numerics. With increasing xq, the L -2 / 3 and L~ 1 / 2 -corrections approach each other. 
The numerical values have been obtained by monitoring the average magnetization 
M defined in 134|) and determining the plateau height, where M jumps from a finite 
value to zero. The slight modulation of the red line arises from the finite numerical 
resolution of this procedure. 

3.6. Numerical results 

To test the analytical predictions derived in the preceding subsections, we have carried 
out a detailed numerical study of the dependence of A and Wq on ko, L and 7. In 
figure [H we show two examples for the dependence of A on the plateau width ko- The 
prediction of the asymptotic maximum principle (jl]) reproduces the qualitative behavior 
of the numerical data but significantly overestimates the value of A. The L~ 2 / 3 finite 
size correction (|23|) derived in the harmonic approximation improves the comparison, 
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Figure 3. Illustration of the power of the sequence length in the correction term for 
fixed relative and absolute plateau width. AA is the numerical value for the growth 
rate A num less the value obtained from the maximum principle, eq. (j4|). For fixed 
relative plateau width, as well as for fixed absolute width, the numerics show the same 
exponent of the sequence length as the corresponding analytical result. 



f(k) 



L-ki L 



k 



Figure 4. Illustration of a fitness landscape with two plateaus. This type of landscape 
is used to investigate the influence of height and (relative) broadness of the plateaus 
on the population fitness A. 



but quantitative agreement is achieved only using the refined expression fl3Tj) . which is 
proportional to L~ x l 2 . 

Figure [2] shows a similar comparison for the critical plateau height Wq. Here the 
prediction (1261) of GH is also included and seen to match the numerical outcome only 
poorly, whereas the MP result with the finite size correction (1311) produces excellent 
agreement. Finally, in left panel of figure [3] we verify that the finite size correction A A 
indeed varies as L~ 2 ^ 3 when L is increased at fixed relative plateau width xq. The right 
panel shows the corresponding L -1 / 2 dependence for fixed absolute plateau width k . 
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4. Fitness landscapes with competing peaks 

4-1. The selection transition 

Since we have validated the analytical results of section [3] via numerical studies, we 
can now apply the analytical theory to the phenomenon of the survival of the flattest 
as explained in the introduction. To be specific, we want to find out whether a broad 
plateau outcompetes a smaller but higher one even in the limit of long sequences. 

In the literature this question has already been discussed to some extent by Schuster 
and Swetina [7J. The question can be answered by investigating a fitness landscape 
consisting of two fitness plateaus at the opposing ends of the Hamming space (see figure 
H]). As was shown in [7J, for small /1 the interference between the two plateaus is negligible 
when they are separated by a few mutational distances. The stationary state of the 
system is therefore to a very good approximation determined by the comparison between 
the population growth rates associated with each of the two plateaus in isolation. 

Observing the center of mass of the population as function of the mutation rate 
and for fixed sequence length, we find two types of transitions. The first one is a jump 
of the population from the higher to the broader plateau, which we will refer to as the 
selection transition [23] taking place at a mutation rate /j, s . The second transition is the 
well-known error threshold taking place at /it r , where the population becomes uniformly 
spread in sequence space. To analyze these transitions, an order parameter is needed. 
A convenient quantity is the population averaged "magnetization" defined by 



If the whole population consists only of master sequences, the magnetization is M = 1. 
If only the inverse master sequence is present, the magnetization becomes M — — 1. For 
a uniform distribution in sequence space (delocalized population) the magnetization is 
M = 0. Thus we can distinguish the qualitatively different states of the population in 
the two plateau landscape by considering the population averaged magnetization M as 
a function of \i. 

As can be seen from figure the selection transition (the jump between the two 
plateaus) is sharp even for finite sequence lengths, whereas the error threshold is a 
continuous transition for finite sequence length and only becomes sharp in the limit of 
infinite sequence length [23]. With growing sequence length the two critical mutation 
rates /x s and /j,t r become smaller and also approach each other until, at a critical sequence 
length L*, the selection transition completely disappears. For sequences longer than 
L*, the population delocalizes directly from the high, narrow peak and the low, broad 
plateau is never substantially populated. 

With the help of the maximum principle, this surprising behavior can be easily 
understood. Using the selection threshold is obtained by equating the population 




(34) 



k=0 
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mean fitness for the two competing peaks, which yields 

Wq — W\ Wq — W\ 




(35) 



2 L/hL -k\- ySk L - kfj ^{Vh-Vh) 



where ko, fci « I has been assumed in the last step. On the other hand, the error 
threshold jJ^ associated with plateau i = 0, 1 is determined by the vanishing of the 
corresponding principal eigenvalue A,, which gives 



The different scaling of the two types of thresholds with sequence length implies that 
for large L the error threshold of the higher peak is encountered before the selection 
threshold, which therefore is no longer observable. The critical sequence length L* 
where the selection transition vanishes can be estimated by equating the approximate 
expressions fl35l) and fl36|) . which yields 



Following [7J |23], in our numerical work we have considered short plateaus, k Q = 1 
and k\ = 2, with relative fitness values u>i/u>o = 0.9, for which (1371) give L* ps 106. 
Comparison with the numerical values for the selection and error thresholds in figure [6] 
shows that this significantly underestimates the value of L*; moreover, the agreement 
between theory and numerics is not substantially improved by using the full expressions 
for the principal eigenvalues Ao and Ai, including the L _1 / 2 -correction derived in 
subsection 13.51 This is not surprising, as the continuum approach developed in section 
[3] cannot be expected to be quantitatively accurate for plateaus sizes of order unity. 

For completeness we mention that for plateau widths scaling with the sequence 
length (such that Xq = k /L and X\ = k±/L are kept fixed as L —>■ oo) the selection 
transition is maintained at a fixed value of 7 [18J. 

4-2. The ancestral distribution 

In addition to the equilibrium population distribution P£ attained at long times, we can 
also consider the ancestral distribution, the equilibrium distribution of the backward time 
process, as introduced by Baake and collaborators [21] • The ancestral distribution 
0^ gives information on the origin of the equilibrium population and is obtained as the 
product of the right eigenvector P£ and the left eigenvector P£* of the mutation-selection 
matrix defined through ([2D, a k ~ P fc * • P£*. 

For the fitness landscape with two competing peaks, we find an ancestral population 
that is either located on one of the plateaus or uniformly distributed in sequence space 
(figure [7J . The transitions beween these states are all of first order. The continuous 
character of the error threshold transition of the equilibrium population distribution, as 
opposed to the discontinuous transition of the ancestral distribution, can be explained 
by the growing mutational pressure affecting the population on the plateau and driving 




L 



(36) 




(37) 
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L = 500 



oTcJos ate a~ui5 ate 01525 ^ 

IMr 

Figure 5. Order parameter M as function of the mutation rate /i per site for a 
fitness landscape with a high plateau at k = and a broad plateau at fc = L. For 
short sequence lengths one can observe a hopping of the population from the higher 
to the broader plateau and then a derealization (left picture). For long sequences, 
one only observes the derealization transition from the higher fitness plateau (right 
picture). The hopping between the plateaus we call the selection transition. It 
takes place at mutation rate The derealization transition, also called error 

threshold, takes place at a mutation rate p,t r - The underlying fitness landscape is 
w fe = 10 • 6(1 - k) + 9 • Q(k -(L- 2)). 





Figure 6. Critical mutation rate [i s of the selection transition and [itr of the error 
threshold for the fitness landscape Wk = 10- 9(1 — k) +9- 6(fc— (L — 2)) as functions of 
the sequence length. The two lines cross at a critical sequence length L* . The selection 
transition is observed only for sequence lengths smaller than L*. The numerical 
data are compared to analytic predictions based on the MP including the L -1 / 2 - 
correction term. As before, the numerical values have been obtained by calculating 
the magnetization of the population and determining for each sequence length the 
mutation rate where the magnetization jumps. 
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equilibrium distribution 



ancestral distribution 




Figure 7. The dominant entries of the equilibrium and ancestral population 
distribution as a function of the mutation rate per site, calculated numerically. The 
underlying fitness landscape is the same two-plateau landscape used in figures [5] and [6] 
The sequence length is chosen as L = 50. Occupation fractions are plotted only for the 
most populated Hamming classes. The two distributions undergo phase transitions at 
the same mutation rates, but at the error threshold the ancestral distribution undergoes 
a discontinuous transition, while for the equilibrium distribution the transition is 
continuous. 



it towards the middle of the Hamming space. Mutations cause the population to "leak 
out" from the plateau. Nevertheless, the individuals maintaining the population and 
compensating for the mutational loss are the ones with highest fitness, which are located 
on the plateau and make up the ancestral distribution. 

Before closing the analysis of plateau-shaped fitness landscapes, we want to mention 
the connection between our description and the popular language of Ising chains or semi- 
infinite Ising models [271 EH]. in the Ising picture, the ancestral distribution becomes the 
bulk distribution on a semi-infinite two-dimensional (spatial or spatio-temporal) lattice, 
and the equilibrium distribution becomes the distribution in the surface layer. This 
analogy can be seen very clearly in the paper by Tarazona [23J, where the different orders 
of the transitions in the two distributions are explained in terms of surface wetting. 

5. Epistasis and the error threshold 

So far, we have discussed robustness of phenotypes using plateau-shaped fitness 
landscapes, which are a special case of the class of epistatic fitness functions. We now 
want to discuss the latter in a more general framework. Epistasis describes the non- 
linear dependence of the fitness function on the number of mismatches k [13J. Every 
additional mismatch is penalized harder (synergistic epistasis of deleterious mutations) 
or less hard (diminishing epistasis) than the previous one. Here we address the effect of 
epistasis on the existence of an error threshold, which is defined for our purposes as a 
singularity in the dependence of the population mean fitness on the mutation rate. In 
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general (but not always, see below) such a singularity is associated with a discontinuous 
jump in the location of the most populated genotype. 

Following Wiehe [H] we consider the class of permutation invariant (Malthusian) 
epistatic fitness functions 

Wk = wq — bk a , (38) 

where k is again the Hamming distance to the master sequence and b > 0. The epistasis 
exponent a takes the value a = 1 in the non-epistatic case, while a > 1 and a < 1 
produces landscapes with synergistic and diminishing epistasis, respectively. For a — > 
( 138]) reduces to the sharp peak landscape = w — b(l — t^o)- It is well known that an 
error threshold exists for a — > 0, but not for a = 1 [6]. Neglecting backward mutations, 
Wiehe argued in [H] that an error threshold emerges whenever a < 1. In the following 
we show that, based on the maximum principle (j3J), the critical value of the epistasis 
exponent below which an error threshold develops is in fact a = 1/2. 

As before, we work in the scaling limit L — > oo and /j, — ^ with the mutation rate 
per sequence 7 = fj,L = const. In order to cast (1381) into the form (j3j) required for the 
application of the maximum principle, we write 

w k = f{%) = ujq — bx a , with x = k/L, b = bL a (39) 

and the limit L — > 00 should be combined with 6 — ^ 0, such that b = const. Since b 
can be interpreted as a kind of selection coefficient, we are thus considering a situation 
where both the mutation rate (per site) and the selection forces are small. Applying the 
maximum principle (HI) to this landscape, the mean fitness A of the population in the 
equilibrium state is given by 

A = max{u> — bx a — 7II — 2\/x(l — x)]\ = max X(x), (40) 

xg[0,l] x£[0,l] 

where X(x) is the function inside the curly brackets. 

To find the condition under which the maximum is attained inside the interval 
x G (0, 1), we set dX/dx = 0, yielding the condition 

1(1 -2x) =ax a ~ 1/2 VT zr x~. (41) 
b 

For a > 1/2 the right hand side is a convex function which vanishes at x = 0,1, with 
an infinite slope at x = 1. As a consequence, there exists always a unique solution 
x c G (0,1) for any value of 7/6, which describes the location of the population for 
L — > 00. The location varies smoothly from x c = for 7/6 — >• to x c — > 1/2 for 
7/6 — > 00, and there is no error threshold. However, for a < 1/2 the right hand side 
diverges at x — 0, and there is no solution for small 7/6. The function X(x) is then 
monotonically decreasing, which implies that the maximum in (1401) is located at the 
boundary point x = over a finite interval of 7/6. Increasing 7/6 the function X(x) 
develops a local maximum, which eventually exceeds the boundary value A(0) = wo — 7. 
At this point the population discontinuously delocalizes to an interior point x c G (0, 1). 
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Figure 8. Magnetization as a function of mutation rate for the fitness landscape 
(|38|) with epistasis exponent a = 0.4 and a — 0.52, respectively. For a — 0.4 the 
magnetization undergoes a discontinuous jump, whereas for a = 0.52, it changes 
smoothly. Calculations have been done for a sequence length of L — 500. 



The error threshold condition is of the form 7/6 = g(a) where the function g(a) is not 
explicitly available. This translates into the expression 

Utr = -j~ = —j— = g{a)bL , (42) 



for the critical mutation rate jj tr . This scaling of /i tr with L was also obtained in 
In the sharp peak limit a — > the threshold occurs at 7/6 = 7/6 = 1, which implies 
that g(0) = 1. On the other hand, for a = 1/2 the expansion of A(x) near x = reads 

X(x) &w -(b- 2 7 )a; 1/2 - 7 x 3/2 , (43) 

which shows that g(l/2) = 1/2. For 7/6 > 1/2 an interior maximum appears at 
x c = (2 — b/^f)/3, which moves continuously away from x = 0. In the language of phase 
transitions, a = 1/2 can thus be viewed as a critical endpoint terminating the line of 
discontinuous phase transitions that occur for a < 1/2. 

These predictions are fully confirmed by numerical calculations for finite sequence 
length. Figure M illustrates the existence of an error threshold for a < 1/2 and its 
absence for a > 1/2 by showing the behavior of the magnetization M as a function of 
7 for two different cases. The magnetization displays a non-analytic jump for a < 1/2 
and varies smoothly for a > 1/2. In figure [9] we show the error threshold as a function 
7/6 = g(a), which interpolates between the limits g(0) = 1 and g(l/2) = 1/2. 

It should have become clear that the special role of a = 1/2 derives from the fact 
that for this value the leading order behavior of the fitness function for small x matches 
that of the "entropic" term ~ — x) in the maximum principle (HI). Since a similar 

term appears also for general alphabet sizes [see fl50l) ]. the considerations of this section 
hold in that case as well. 



Robustness and epistasis in mutation- selection models 17 

7 




a 



0.1 0.2 0.3 0.4 0.5 



Figure 9. Numerically determined phase diagram for the epistatic fitness landscape 
(I38p . At the thick line the population undergoes a first order phase transition from a 
state localized at x c — (below the line) to a delocalized state x c > (above the line). 
This line terminates in a second order phase transition at a — 1/2. The deviation 
from the prediction 7/6 = g(l/2) = 1/2 at a = 1/2 is due to finite sequence length 
corrections. For all larger values of the epistasis exponent, a > 1/2, the population 
changes smoothly. Calculations have been performed for a sequence lengths of L = 750. 
The slight modulation of the line is due to the finite numerical step size. 

6. Conclusions 

In this paper we discussed the properties of epistatic fitness landscapes, with particular 
emphasis on mesa landscapes describing mutational robustness of phenotypes, which 
have been studied previously in the context of regulatory motifs [31 Hj. As population 
evolution model we used the continuous time Crow-Kimura model, which is a 
quasispecies model for asexual and haploid organisms, and analysed its stationary states 
for sequences consisting of two letters. As explained in Appendix B, it is straightforward 
to generalize our results to sequence alphabets of general size. Similarly, the extension to 
discrete time dynamics can be carried out by replacing the Malthusian fitness landscape 
Wk by its Wrightian counterpart W\~ ~ exp(u> fc ) [SI CEE1 [19] . 

We reviewed two existing approaches [31 [TS1 IE] to this problem and explained 
the discrepancy between their predictions by extending the approach of Gerland and 
Hwa [3] beyond the harmonic approximation. Based on a quantum mechanical analogy 
we derived a novel finite size correction term to the maximum principle of [IB] , which 
significantly improves the agreement with numerical calculations. Our central result is 
that the relative number of tolerable mismatches xq = ko/L is the relevant parameter 
for the fitness effect of mutational robustness, and we provide accurate formulae for its 
quantitative evaluation. As a consequence, we showed that the selection transition first 
described by Schuster and Swetina [7] disappears for long sequences. 

Finally, in section we discussed more general forms of epistatic fitness landscapes 
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with regard to the existence of an error threshold. Based on the results of [THl ETJ we 
improved on earlier work [13] and showed that diminishing epistasis [a < 1 in the fitness 
function (1581 ] is not a sufficient condition for an error threshold to occur. 
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Appendix A: The large deviations approach 

We start by symmetrizing the eigenvalue problem (j2J). The discrete analogue of the 
transformation (fl4j) is 

L X 1/2 



Qk=[ k ) p k, (44) 



which leads to 



A<2fc = (to* - i)Q k + - k)(k + 1) Q k+1 + fiy/(L - k + VjkQk-i. (45) 

Following [20], we now perform the continuum limit by making a large deviations ansatz 
for Q k , 

Qk = QxL = ip{x) = exp[-e _1 ?i(x)] (46) 
with e = l/L. Inserting this into (1451) one finds 



(A - f{x) + 7 )V> = 2 7v /x(l - x) cosh.[u'}ip . (47) 

Cancelling ip on both sides yields a Hamilton- Jacobi equation for the 'action' u(x), 
with u' = du/dx playing the role of a canonical momentum [20] . In order to cast ( 1471) 
into the form of a Schrodinger equation, we expand the momentum-dependent factor to 
quadratic order, cosh(it') ~ 1 + (w') 2 /2, and make use of the relation 

M , = e ^~5?' (48) 

which follows from (|46p to leading order in e. Inserting this into (1471) results in 



Appendix B: General alphabet size 

Here we show how our results generalize to the case where the symbols in the genetic 
sequence are taken from an alphabet of A > 2 letters (for nucleotide sequences A = 4). 
We assume a uniform point mutation rate fi connecting any two of the A possible states 
of a site in the sequence, and a fitness function Wj~ that depends on the relative number of 
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mismatches according to <^j. It is then straightforward to see that the basic eigenvalue 
problem (j2J) generalizes to 



(a - w k )p; 

= (A- 1) 7 



(49) 



^±ip; +1 + (L - k + 1)P* 



-i 



L — k 



k 



A-l 



TD* 



Applying the results of [HI ESI 130] to this problem we find that, asymptotically for large 
L, the principal eigenvalue is given by the maximum principle 



A = max < fix) - (A - 1)7 

xe[o,i] 



_ (A-2)x \ _ 2y/x(l-x) 
A-l J y/A^l 



(50) 



For the case of the mesa landscape flU]) this implies that the population is localized near 
the optimal genotype for w > Wq with 



Wn 



j(A-l) 



_ (A - 2)x _ 2y/x (l-x ) 
A-l ^A^l 



(51) 



The right hand side is a monotonically decreasing function of xq which vanishes at 
Xq = 1 — 1/A. For fixed Xq it is an increasing function of A. 
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